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We analyze deficiencies of commonly used Coulomb approximations in Generalized Born solvation 
energy calculation models and report a development of a new fast surface-based method (FSBE) for 
numerical calculations of the solvation energy of biomolecules with charged groups. The procedure 
is only a few percents wrong for molecular configurations of arbitrary sizes, provides explicit values 
for the reaction field potential at any point of the molecular interior, water polarization at the 
surface of the molecule, both the solvation energy value and its derivatives suitable for Molecular 
Dynamics (MD) simulations. The method works well both for large and small molecules and thus 
gives stable energy differences for quantities such as solvation energies contributions to a molecular 
complex formation. 



I. INTRODUCTION. 

Solvent plays an essential role in biophysics in de- 
termining the electrostatic potential energy of proteins, 
small molecules and protein-ligand complexes. Solva- 
tion energy is a major contribution to the protein fold- 
ing problem and to ligand binding energy calculations. 
In the latter case it is the interaction, which is pretty 
much responsible for binding selectivity [14j [28]. Large 
scale Molecular Dynamics (MD) simulations [TJ [18j [23] 
or industrial-scale calculations of the solvation energy 
in drug discovery applications require a fast method ca- 
pable of dealing with arbitrary molecular geometries of 
molecules of vastly different sizes within a single, fast, 
numerically robust framework. 

A solvation energy calculation for a molecule-sized ob- 
ject has always been and still is a challenging problem. 
The most accurate approach is, apparently, a large scale 
MD simulation [24j [25] of the body of interest immersed 
in a tank of water molecules in a realistic force field or 
even within quantum mechanical settings. Although be- 
ing ideologically correct such calculations are time con- 
suming and pose a number of specific problems stem- 
ming, e.g. from long relaxation times of water clusters. 
One possible way to bridge such "simulation gap" is to 
employ different types of continuous solvation models. 
Fortunately, water is characterized by a very large value 
of dielectric constant and therefore to a large extent the 
reaction field of water molecules has a collective nature. 
Although realistic properties of molecular interactions 
depend both on short-scale water molecules alignment 
and on their long-range dipole-dipole interactions at the 
same time [7j [8], purely electrostatic models, such as 
Poisson-Boltzmann equation solvers [2j [29], turned out 
to be very successful in various applications. 

Even within the realm of continuous electrostatic mod- 
els there are numerous approaches in use to calculate the 
electrostatic contribution to solvation energies. Popu- 
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lar techniques span from finite element methods (FEM, 
0[l[5l[6l[Il[2l[35l[37])to multiple variations of Gener- 
alized Born (GB) approximations [TT1 [Tol [20l [22l [26l f28l 
E01 EU [34]. A numerical FEM solution to the Poisson- 
Boltzmann equation (PBE) is a formally fast (the cal- 
culation time and memory scale ex TV, with N being the 
number of particles in the system) and is a rigorous at- 
tempt to solve the electrostatics problem. On the other 
hand GB approximations are practically fast, in spite 
of the fact that it normally takes 0(N 2 ) operations to 
calculate GB energy. Unfortunately GB approximations 
are very rough and that is why GB calculations work well 
only for small and medium sized molecules, whereas FEM 
methods can, although at expense of numerical complex- 
ity, be applied to very large systems. The particular 
boundary between the applicability of the two methods 
is vague and depends, in terms of speed, on the details 
of the methods realization, and, in terms of accuracy, on 
the system geometry (see below). 

In this Paper we report a development of a new fast 
surface-based method (FSBE) for numerical calculations 
of the solvation energy of biomolecules with charged 
groups. First we elucidate physical nature of commonly 
used GB models, identify the variational principle behind 
and discharge the so called Coulomb approximation. As a 
result we suggest a new computational procedure, which 
is only a few percents wrong for any molecular config- 
urations of arbitrary sizes, gives explicit values for the 
reaction field potential at any point inside a molecule, 
characterizes the water polarization charge density on the 
molecule interfaces. The approach reported here is suit- 
able both for the solvation energy and its derivatives cal- 
culation for Molecular Dynamics (MD) simulations. The 
method works well both for large and small molecules and 
thus gives stable energy differences for quantities such as 
solvation energies of molecular complexes formation. 

An important side effect of our studies is a comparative 
research of various GB approximations. We distinguish 
between the volume and surface based approaches to cal- 
culate the Born radii of the charges and demonstrate 
that only the latter can be trusted. The reason is that 
any practical way of volume overlaps integrals calcula- 
tion implies some sort of weak atomic overlap approxima- 



tion and hence effectively leaves many unphysically small 
water-filled cavities within the molecules. This leads to 
an overestimation of both the molecular volume and the 
dielectric constant of the molecular interior. Both factors 
essentially disrupt accurate descreening calculations and 
often lead to completely unrealistic results. 

The manuscript is organized as follows. Section II pro- 
vides an overview of continuous solvation energy calcu- 
lation methods. We compare exact Poisson-Boltzmann 
equation solvers to approximate GB models and highlight 
deficiencies of commonly used Coulomb approximations. 
In the subsequent Section III we represent the idea of a 
new fast molecular surface based method and estimate 
its accuracy for a number of exactly solvable cases. In 
Section IV we discuss important numerical implementa- 
tion details and, at last, in Section V, we demonstrate 
performance of the method in realistic modeling exam- 
ples. At the end of the presentation we hint how explicit 
expressions for the surface charge densities and the reac- 
tion field potential may help building O(N) methods for 
approximate solvation energies calculations [TO] . 
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Figure 1: Solvation energy calculation problem setup: 
schematic representation of a macromolecule (see the expla- 
nations in the text). 



II. OVERVIEW OF EXISTING METHODS. 

To elucidate the nature of approximations and limita- 
tions of GB family of approaches it is instructive to start 
from the basics physics. To find the polar contribution 
to the solvation energy in a continuous solvation model, 
Esj one should solve the Poisson equation 

A<p(r) = -47rp(r) (1) 

for the potential y?(r) generated by the charge density 

p(r) = 5><5(r-r,), (2) 

i 

defined by the atoms placed at the positions r^, and the 
boundary conditions at the molecules surfaces and spatial 
infinity. 

There are various ways to calculate the potential (p(r). 
The most practical approach is to use some sort of finite 
elements method (FEM), which can be both in volume 
and boundary grids incarnations (see e.g. [4j |5j [6j [16j 
HE El OH GEl [37]). The boundary grid based methods 
are often more practical and aside from subtle details 
are equivalent to Surface Electrostatic Solvation (SES) 
models. A typical SES-water model can be considered 
as an alternative to discretization of the volume and is 
given by the solution of the following integral equation 

2na 3 (r) + f dfaj (r') = -q 3 (3) 

Jr w |r r | |r Tj\ 

for the polarization charges surface density aj (r) at the 
point r on the molecule's surface induced by the protein 
charges q 3 - as shown on Fig. [T] Here df is the element of 
the molecular surface at a point r', n is the unit normal 



to the surface at the point r. The exact formula for 
solvation energy is then: 

(E S ) ex = \Y,1^ r ^ ( 4 ) 

i 

where 

stands for the so called reaction field potential, produced 
by the water polarization charges on the boundary of the 
molecule Tw The total electric potential consists of the 
two parts: 

(p(r) = <p (r) + y>i(r), (6) 

where 

3 = 1 1 Jl 

is the potential of the charges in vacuum, i.e. in the ab- 
sence of the water molecules. Since water is characterized 
by a large value of the dielectric constant, ew ~ 80 ^> 1, 
to a good accuracy the electric potential vanishes inside 
the water bulk so that 

<^(r)|r w =0 (8) 

on the boundaries. The model implies that the dielectric 
constant of the liquid is infinitely large, whereas the di- 
electric constant of the molecules interior is 1. Although 
the method is fairly easy to implement, it is also not 
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very practical: in realistic applications involving large 
molecules the calculation is memory consuming, slow and 
not very stable with respect to small changes in the sur- 
face elements positions and orientations. The latter cir- 
cumstance also means that both FEM and SES methods 
often fail to provide smooth derivatives of the solvation 
energies suitable for MD studies of bio-molecules. 

A very well known alternative to solving the Poisson- 
Boltzmann equation directly is to use generalized Born 
(GB) approximation, which is a fast, simple, qualitatively 
correct and numerically stable method for macromolecu- 
lar solvation effects calculations Hfl E2J El [34] . The 
method is based on the following ad hoc. approximate 
expression for the full electrostatic energy E e \ for sys- 
tem of charges charges qi located within the surface Yw 
separating the molecule from the water environment: 



(9) 



spite of being only a very rough approximation, GB mod- 
els are widely used in practical simulations. 

Common deficiencies of GB approximation are very 
well known. Consider, e.g., a single charge q fixed at a 
distance r from the center of a spherical molecule of a 
radius a. Eqs. (10)-(12) immediately yield: 
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(13) 



(14) 



On the other hand the problem is simple and can be 
solved exactly both for the reaction field potential [T7| 

[32]: 



(15) 



The notations used in the expression are illustrated on 
Figjl] Here the indices z,j = 1,...,7V enumerate the 
charges, N is the number of charges, = r ( 
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r^ |, Yi is the radius- vector of a charge i (i-th 
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(10) 



ep and ew are dielectric constants for within the 
molecule interiors and water, correspondingly. The fac- 
tor fcB (Uj) is commonly (although not always) defined 
as 



fen (nj) = \r% + R B iRBjexp (-r^/4,R Bi R Bj )] 



1/2 



(11) 

The effective Born radii Rb% of the ions are calculated 
according to 



4tt J w sj 4tt J rv 



(12) 



J_ = J_ f lM = ±f ^ 
Rb% 

where Si = |s^|, = r' — r^. In its volume integral 

representation Eq. ( 12 ) assumes the integration over the 



water bulk W . which can be easily transformed to an 
equivalent boundary integral form in a standard way with 
the help of the Gauss theorem [T3] . 

Various models are used to define molecular surfaces 
and volumes. Normally the molecule volume is approxi- 
mated as a set of spheres of specified radii a^, the individ- 
ual ions Born radii, centered at the points of the charges 
locations and thus characterizing water cavities associ- 
ated with the ions in the solute. Therefore a complete 
GB model should also include a set of fitting parameters 
di. The specific values of the model radii are either set 
to the atomic van der Waals radii, or (better) trained to 
reproduce experimental values of the polar part of small 
molecules solvation energies whenever it is possible. In 



(f j = Tj/rj), and the solvation energy 
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(16) 



a 2 — 2r ? -i\- 



for an arbitrary number of the charges within the sphere. 
The solution has been long advocated by Kirkwood [19j 
33jand takes especially simple form for a single charge 
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(17) 



The approximate GB solution (fl4| fails to reproduce the 



exact result (17) for the solvation energy of an ion within 
a spherical cavity, as shown on Fig(2] The solvation en- 
ergy and hence the Born radius are in a good agreement 
with the exact result if the charge is close to the cavity 
center and are off by a factor of 2 if the charge is next to 
the molecular surface. Realistic biomacromolecules are 
large and most of their charges are close to molecular 
surfaces. In the very same time the GB approximation 
in its most commonly accepted form fails exactly next to 
the molecular surface. This means that there is no way 
to "train" the seed values of the Born radii to reproduce 
the solvation energies of both small and large molecules. 
This also means that GB models in the standard form 
can not predict well solvation energy contributions to lig- 
and binding free energies in drug discovery applications. 
Indeed, drug binding affinity depends on the solvation 
energy difference between a protein-ligand complex and 
the protein-ligand pair separated at infinity. The proteins 
are large and ligands are normally small molecules with 
all the substantial charges of the protein-ligand complex 
arranged close to the (large) protein surface. 

The reason why GB approaches fail becomes clear from 
comparison with the exact expression 



(E s ) ex = (E S ) GB + AE S < (Es) 

GB i 



(18) 
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Figure 2: Ratio of GB solvation energy to exact one for the 
model spherical "protein" of a radius a as a function of the 
charge position r from the center of the sphere. 



where 



p 07T 



<o, 



where the integration is performed over the molecule in- 
terior P. GB approximation accounts for the electro- 
static energy of the polarization charges (reaction field) 
incorrectly and, in fact, overestimates it. Eq. ( [l8] ) sug- 
gests that GB is nothing else but a variational calcula- 
tion of the solvation energy. The "probe function" (10) is 



widely tested and trusted, whereas specific recipes for the 
Born radii calculations can still be different. The popular 
choice, Eq. ( [l2| ), corresponds to the so called Coulomb 
approximation (CA, see [3] and the refs. therein for a 
review). CA does not follow from any first principles 
and puts severe limitation on applications of GB models. 
Up to date there have been a few sound attempts to go 
beyond CA and obtain better recipes for the Born radii 
as discussed, e.g., in [13j [27] . In what follows we dwell 
into the physics behind the Born radii calculations and 
generate a whole family of approximations for molecular 
electrostatics. 



III. HOW TO FIND BORN RADII? 

In this section we part from C A and demonstrate a new 
way to calculate the polar part of the solvation energy. 
The practical goal is to combine the accuracy of FEM or 
SES models with the speed and numerical stability of GB 
approximation. To prove this is possible we identify GB 
solution as a possible variational solution of the Poisson 
equation ([!]) . Given a set of known positions of the atom 
charges, we suggest the following GB-like anzatz for the 
reaction field potential tp\\ 



<Pi(r) 



(r-Vjf + R^Rj 



(19) 



wherei?(r) is the variational function, Rj = R(rj). The 
true solution of the electrostatics problem provides the 



minimum to the functional: 



G 2 [R(v)} = J p dV^(V Vi r. 

Since the potential vanishes at the molecule boundary 
Eq. (|8| suggests a very simple boundary condition for 
the variational function R(r): R(r) |p w = 0. 

The potential in the form of Eq. dl9| is an approxi- 



mation already. The best possible function R (r) should 
provide the minimum to the functional G2. To find such 
a solution may be an interesting problem in itself. Never- 
theless it is not practically: optimization of the functional 
G2 is roughly as easy (or difficult) as to find the exact 
solution of the Poisson equation. To avoid this unneces- 
sary procedure of the functional minimization we suggest 
instead a specific form of the function R (r) 
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[R(t) 
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~d 3 r', 



(20) 



in the "classic" volume integration form, or, equivalently, 
in the surface integration form 



- = -f 



(21) 



for each of the charges. Here si = |s^|, = r' — r^, and 
the polar part of the solvation energy (the reaction field 
energy) is given by a Kirkwood like expression 
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(22) 



with = f(ri,rj) = 



4 



■R(ri)R( rj ) = 



RjRj 



Although at a first glance FSBE approach does not 
seem to be very different from GB approximation, the 
solution (19) is a much better approximation to the solu- 



tion of the original electrostatic problem. To see that let 
us turn back to the example of a charge confined within a 
spherical cavity of radius a. The new improved Eq. (20) 
for the "generalized" Born radius gives 



R (r) = (a 2 



/a, 



(23) 



which, after inserting into Eq. ( fl9| gives the exact results 
for the reaction field potential ( |15| ) and the solvation en- 
ergy of the point charge ( [l6] ) within the sphere. It can 
be further shown that FSBE approach is exact for arbi- 
trary configuration of charges confined within a spherical 
cavity of arbitrary size. This means FBSE is exact both 
for ions next to a large protein boundary and in a center 
of a small sphere representing a single ion. The FSBE 
gives also the exact result for arbitrary configuration of 
multiple charges next to the spherical water cavern inside 
a large protein. 

Our direct interpretation of the reaction field potential 
helps us to find the polarization surface charge density as 
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Figure 3: Ratio of FSBE solvation energy to exact value for 
one charge inside protein in the form of a layer with thickness 
L (the lower curve). The upper curve describe the result of 
the improved approach FSBEi (see below). 



at the interface boundary. Indeed, the standard form of 
the electrostatics boundary condition for the electrostatic 
potential reads: 



1 dip 
4-7T dn ' 



where 
<p(r') 



<y?0 + <Pl 



2* 

3 



/(r',r,) 



is the full electrostatic potential. Next to the boundary 
(r' — > Tw) R (r') w 2h — > 0, where h is the distance from 
a given point to the surface. Combining the expressions 
above we obtain: 



|3 ' 



(24) 




Figure 4: Ratio of FSBE solvation energy to exact value for 
one charge inside the corner between two perpendicular infi- 
nite walls (the lower curve). The upper curve describes the 
result of the improved approach FSBEI (see below). 
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where z = z/L. Once again, to characterize the difference 
between the approximate FSBE and the exact results we 
plotted the ratio of (Es)fsbe ^° ^ ne exac ^ solvation en- 
ergy (Es) ex on Fig|3j As in our spherical cavity example 
above the two results coincide at the dielectric boundary 
(as it should be) and deviate from each other in the cen- 
ter of the layer. The discrepancy does not exceed 9%, 
which is nothing compared with the factor of 2 in the 
case of the standard GB approximation. 

Another challenging case is the calculation for a single 
charge q placed within a corner made of two perpendicu- 
lar infinite walls (the u xz" and "yz" planes). Once again, 
our FSBE result 



Note, that the standard GB approach may, in principle, 
also be used to calculate as- Nevertheless such an ap- 
proximation would not be good since GB approximation 
for R(r) is twice as small than that of the exact result 

FSBE can not, of course, be exact for an arbitrary 



molecule geometry. Eqs. (20) and (22) are certainly only 
approximate. To see the limitations of the approach we 
explored various exactly solvable charges configurations. 
Consider the first example: a plain layer-like "molecule" 
(or membrane) of the thickness L surrounded by the con- 
tinuous water on both sides with a charge q placed inside 
the layer at the distance z from one of the water interface 
planes. The exact result for solvation energy is [T7J [32] 



(Es) ex = q- 



f 

Jo 



dk 



sinh (kz) sinh (k (L — z)) 1 
sinh (kL) ~ 2 



(25) 

Eqs. ([20]) and \22\ be used to find FSBE approximation 



for the solvation energy 



FSBE 



yjl — | (sin (p cos (p) z 
4r sin cp cos cp 



where cp is the azimuthal angle between the position of 
a charge and the "xz" plane, r is the distance from the 
charge and "z" axes (the intersection of the walls). The 
result should be compared with the exact solvation en- 
ergy 



(Es) e 



-q 



2 sm <p + cos if — sm cp cos ip 
4r sin cp cos ip 



Once again, the ratio of the two energies is plotted on 
Figj4| The difference is no more than 6% in the center of 
the system and disappears at the corner boundaries (as 
it should be). 

The presented results prove that Eqs. (20) and ([22]) 



defining FSBE approximation do provide a fairly good 
solution of the electrostatic problem in various geome- 
tries. Whenever a charge is placed close to an interface 
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Position of an atom, (A, counted from the protein center) 



Figure 5: Born radii calculated for an ion placed at different 
positions inside a model "protein" made of 1000 carbon atoms. 



boundary, FSBE becomes exact; for charges placed at 
the central regions of a large molecule the error is about 
10%, which is fair and often not very important, since 
most of the charges in biomolecules are located in a layer 
next to molecular surfaces. This error can be lowered up 
to 2% by further variational improvements of FSBE (see 
below). 

Before we proceed to explicit description of the method 
implementation, let us take a note on volume and surface 
integrals methods for Born radii calculations. Practical 
applications of Generalized Born models are further com- 
plicated by various approximations introduced for vol- 
ume (or surface) integrals calculations. Since direct cal- 
culations are often prohibitively time consuming, the in- 
tegrals are often estimated in various sort of pair approx- 
imations with subsequent removal of the atom overlaps 
etc. Obviously atoms in biomolecules are fairly densely 
packed and the approximation lead to wrong molecular 
volumes and very wrong (even negative(!)) values for the 
Born Radii for every atom. 

Physically speaking Born radii quantitatively show a 
degree to which an atom is "buried" within a molecule, 
such as a protein. Figj5] gives a simple idea to which 
extent GB can even be used for description of solvation 
energies of a simple, model spherical "protein" molecule 
built of approx. 1000 carbon atoms. The red squares 
represent Born Radii as a function of an ion position off 
the center of the "protein". The values were obtained 
using our own implementation of AGBNP method |12| , 
one of the best realizations of GB procedures available in 
the literature. The yellow curve represents exact result 
for a spherical protein. As one can see, AGBNP results 
fail grow enough inwards and saturates at a very small 
value at the "protein" center. 

The explanation is the following: AGBNP (and for 
that reason practically any other GB model based on 
volume integrals approximations) implies a certain im- 
plicit approximation for the shape of molecular surface. 



Since the model equations employed for the atomic over- 
lap integrals do not provide a direct interpretation, it 
turns out that the overlap integrals are often not exact. 
Physically this means that there are effectively numerous 
water filled cavities of nonphysically small sizes assumed 
inside the protein. The cavities are so small that can 
not hold a single water molecule inside, though represent 
(within the same model) a medium with high dielectric 
constant, effectively increase the dielectric constant of the 
"protein" and therefore decrease the value of the Born 
radii. To check the hypotheses we implemented a sim- 
ple algorithm to search for the water filled cavities and 
remove them (to a certain adjustable extent). The re- 
sult is represented by the blue circles and shows a clear 
improvement towards reproducing the exact analytical 
result. The simple exercise shows that volume integral 
based Born models overestimate the dielectric constant 
within the molecule and may easily lead to a number of 
undesired unphysical issues. In practice any approach 
based on a calculation of surface integrals for Born radii 
has much better chances to yield meaningful results. 

Fig(5] demonstrates another feature of Born approxi- 
mations. As discussed earlier CA fails at the protein 
boundary and gives the Born radius which is twice the 
exact result (see the dots on the right compared to the 
yellow line). This is a genuine problem of CA and can be 
solved by, e.g. switching to FSBE expressions for Born 
radii. 

In principle, Eq. (20) can be used to calculate Born 
radii directly. Unfortunately such a procedure is too 
slow for realistic molecules with typical number of atoms 
N ~ 10 4 . Below we will show that FSBE in the form 
of Eq. (21) yields to a much better GB solvation energy 
calculation implementation. Since the solvation energy 
is often used in MD simulations, we need also analyti- 
cal and easily implementable prescriptions for the forces 
calculations, i.e. energy derivatives with respect to the 
atomic positions: 



dE s 
<9r 7 - 



Qk^jk 



It. 



qiqk 

^ i,k (z^) 



Ri 



dRj 



(26) 



Let us show how our surface integral representation of the 
Born radii (21 ) let us to express the forces in terms of the 



surface integrals. To calculate the derivative dRi/drj we 
shift the atom j with coordinates rj by a small value drj 
and observe how the surface elements df are affected by 
the atom move. Then the molecule volume changes by 
the value dV = dr^d? , which lets us calculate the Born 



radius change using the Eq.(21) as follows 



dRj 



47T 

dRj 
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dRj 



(27) 



(28) 
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Solvation energy of an atom pair 




distance between the atoms, A 



Figure 6: Solvation energy of a diatomic molecule (in units of 
1389 kJ-A/ (mol • e 2 ) ) in frames of FSBE approach. The green 
spheres at the inset represent ions, the blue crosses represent 
the surface points that were used in calculation. 



where T J W represents the part of the molecular surface 
influenced by the atom j. In the following section we 
show how GB implementation defined by Eqs. (21), (26) 
and (27) performs in a few model and realistic situations. 



IV. RESULTS AND DISCUSSION. 

FSBE is not mere another method for quantitatively 
correct molecular modeling calculations. In what follows 
shortly we will show that FSBE calculations have a num- 
ber of important properties besides its speed. To see 
that let us consider a few model calculations to show the 
method performance in a number of simple but challeng- 
ing limiting cases. 

A diatomic molecule is the simplest but the at the 
same time conceptually important example of a realis- 
tic solvation energy calculation. The trick is that any 
reasonable solvation energy model gives exact value for 
a single atom. Depending on the radii of and the dis- 
tances between the atoms the solvation energy of a pair 
may be a very good test of a solvation energy model and 
transferability of its parameters. Fig(6] shows the FSBE 
calculated solvation energies for a pair of model ions with 
similar (red curve) and opposite (blue) charges of 1/2 
atomic units each. The results are pleasing and easy to 
understand. At infinite separation both curves saturate 
at —0.125, which is the correct Born solvation energy 
limit in units of 1389 A: J • A/ (mol • e 2 ) for a pair of the 
charges corresponding to bare radii 2. If the total charge 
is (the blue curve) , at r = we have Es = as it should 
be for a neutral system. If the total charge is 2 x 0.5 = 1 
(the red curve), then at r = we have Es = —0.25, as 
it should be for a combined charge within the sphere of 
radius 2. 

Although the asymptotic values on the graph are fine, 
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Figure 7: Solvation energy of a diatomic molecule (in units of 
1389 /c J- k/(mol-e 2 )) for zero total charge. 
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Figure 8: Solvation energy of a diatomic molecule of total 
charge 1 in units of 1389 k J • A /(mol • e 2 ). 



this does not mean that the whole curve is reproduced 
correctly. To compare our approach with the true solu- 
tion of the electrostatics problem and standard GB mod- 
els we performed the calculation of the diatomic system 
by solving the Poisson equation exactly and with the 
help of by two "classic" GB models (that of HCT and 
AGBNP). The results for a diatomic molecule with zero 
total charge are represented on Fig(7| (charges of ions are 
opposite and equal 1/2 and —1/2 ). 

The electrostatic part of the solvation energy corre- 
sponds to the blue curve of the previous graph and is cal- 
culated either by a (surface-electrostatic) Poisson equa- 
tion solver (blue), FSBE (cyan), AGBNP (yellow) and 
HCT GB model (yellow). As it is clear from here, all 
the approaches give very similar results for the "small" 
molecule and are practically indistinguishable. Indeed, 
it is well known that practically any sort of GB approxi- 
mation gives good results for solvation energies of small 
molecules. 

The difference between FSBE method and "classic" 
GB approaches and its relation to the exact solution be- 
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Figure 9: Solvat ion energy of a cluster of total charge (units 
as in Figs |6|7|8 >. 



comes more obvious if we consider a charged diatomic 
molecule, namely, a molecular ion with total charge, say, 
1 placed on one of the atoms (see Figj8|. The exact (blue) 
and FSBE (cyan), once again, are both in agreement with 
each other, whereas both "classic" GB approaches, HCT 
and AGBNP fail to recover correct asymptotic value at 
zero inter-atomic separation. The latter difference be- 
tween GB solutions and the exact value of the solvation 
energy is not important for small molecules (low atom 
density) but is extremely important for macromolecules 
simulations and ligand binding calculations. 

Binding energy calculations of a small molecule to a 
large protein often pose a difficult problem: a method 
for molecular electrostatic energy calculation should work 
well both for the protein ligand complex, the protein and 
the ligand at infinite separation. The protein and the 
complex are normally large molecules, whereas the lig- 
and is, by definition, small. Not every computational 
approach for the solvation energy calculation is fit for 
the job though. To elucidate the nature of the problems 
at hand we performed another model calculation. First 
we prepared a spherical "protein" of a large (but realis- 
tic) radius. Then we placed a single-atom ligand with a 
charge at a given distance from the "protein" center as 
shown at the insets to Figures [9] and [lOj Then we calcu- 
lated the solvation energy of the system as a function of 
the ligand distance both when the protein is neutral and 
charged (in the latter case the protein charge was taken 
opposite to that of the "ligand") 

Once again we used four different methods for the elec- 
trostatic contribution to the solvation energy calculation: 
a Poisson equation solver (in its surface electrostatic in- 
carnation, blue), FSBE (cyan) and the two "classic" GB 
methods, based on the Coulomb approximation: HCT 
(magenta) and AGBNP (yellow). Fig. [9] corresponds 
to an overall electrically neutral cluster and shows abso- 
lute deficiency of HCT approach deep enough inside the 
"protein". The problem is caused by unrealistic assump- 
tions with regard to the overlap integrals calculations is 
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Figure 10: Solvation energy of a cluster of total charge 1 (units 
as in Figs |6|7|8 >. 



occurs pretty frequently in realistic proteins. AGBNP 
method represents one of the latest and possibly the best 
among GB approaches. In fact the method is specifi- 
cally designed to account for the atoms overlap better 
and ease the problem. However, AGBNP is based on 
Coulomb approximation and thus fails to recover cor- 
rect behavior of the solvation energy close to the "pro- 
tein" boundary: AGBNP energy is off by a large number 
from both FSBE and the exact solution. Remarkably, 
the FSBE and Poisson solutions agree very well every- 
where. Fig|lO] shows the same calculation for a charged 
model "protein-ligand" complex. Once again, HCT fails 
entirely, AGBNP does not work properly at the "protein" 
boundary and both the Poisson solver and FSBE agree 
very well, though FSBE does not require iterations and 
hence is about one order of magnitude faster than a FEM 
Poisson equation solver. 

The results presented in this Section so far may be fine 
but concern only a few oversimplified examples produced 
for model systems with idealized geometries. To judge on 
actual performance of the method we turn to a practically 
interesting realistic system: solvation energy calculations 
for iV8-neuraminidase protein (pdb accession code 2htl). 
The molecule is composed of 387 amino acids and, af- 
ter all the hydrogen atoms added, has 5866 atoms. The 
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results of the calculations are represented on Fig. 
The horizontal axis represents the Born radii obtained 
"exactly" by solving surface boundary condition version 
of the Poisson equation as described by Eq. ([3]). The 
vertical axis shows the Born radii subsequently obtained 
by "standard" CA GB method in its surface incarnation 
rtl2l), FSBE and our in house realization of AGBNP. Both 



the surface Born and FSBE calculations were performed 
using the same surface generated using the same set of 
(realistic) atom radii. The solid dots very next to the di- 
agonal correspond to FSBE results. The values obtained 
with the standard Generalized Born approximation are 
depicted by the turned crosses and generally lay above 
the "exact" results. At last, the AGBNP results are given 
by the crosses at the bottom of the Figure. 
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The results of the calculation support every statement 
we made and hopes we put in designing FSBE method. 
AGBNP does not work well since its pair approxima- 
tion to the overlap integrals estimation does not hold 
for a densely packed atom ensemble, such as a realistic 
protein exactly in the same way as it happened in our 
model spherical "protein" calculation discussed earlier in 
this paper and presented on Fig. |5j It is not a spe- 
cific AGBNP fault, in fact any method based on pairwise 
descreening estimations would perform similarly. FSBE 
appears to fair very well especially when the Born radii 
are small which is indicative to atoms next to the protein 
surface, where normally all the ions are and most of inter- 
esting interactions, such as protein-ligand coupling occur. 
Standard GB in CA fails to reproduce Born radii values 
smaller than 10 A. In fact the radii calculated in CA are 
two times larger than those obtained with FSBE or the 
exact values. It is exactly the behavior we expected from 
our earlier sphere model discussion (see Fig. [2]). The 
Figure also shows that neither FSBE results are perfect. 
Nevertheless FSBE is clearly superior to surface GB in 
CA, provides better both quantitative (at low Born radii 
values) and qualitative agreement with the exact results. 
The apparent deficiency of the method for large Born 
radii is also explainable: large Rb correspond to deeply 
buried atoms, which is exactly the situation when FSBE 
results deviate from the exact solution most. We note 
that FEM such as SES are merely attempts to solve elec- 
trostatics problem in a complicated molecular geometry 
and may be sometimes produce wrong energies due to its 
own method specific problems. 



V. CONCLUSIONS 

The results and the analysis above suggest that our 
FSBE approach represents a fast and fairly accurate ap- 
proximation to the Poisson equation solution. FSBE ap- 
proach does not rely on Coulomb approximation (CA) 
and is shown to work well both for small molecules and 
large molecular clusters involving molecules of very dif- 
ferent sizes. Therefore, FSBE has a potential to compute 
solvation energies with a single transferable set of GB pa- 
rameters capable of describing correct dissociation limit 
of large and small molecules on the same footing. 

FSBE is conceptually simple and shares the best of 
the two words: the calculation speed and smoothness of 
the energy surface of GB models and accuracy of FEM. 
Therefore the approximation should become a weapon 
of choice for a (relatively) fast calculation of solvation 
energies in modeling. FSBE is not a rigorous variational 
solution to the Poisson equation and can therefore be 
further improved. Neither FSBE is the only possible way 
to get rid of CA. As suggested earlier, both FSBE and 
even "classic" GB can be viewed as a variational approach 
with, e.g., single-parameter probe function of the kind: 




+ + 4* + + + :&+ ++ + +++ + 




Figure 11: Born radii calculation as a comparison of different 
GB methods performance vs. "exact" Poisson equation solver. 
The atoms position were taken from a crystallized structure of 
NS neuraminidase (pdb accession code 2htl). The dots, the 
turned and the standard crosses correspond to FSBE, surface 
GB in CA and a volume integral with pair overlaps estimation 
(AGBNP). 



where a is the variational parameter, and C a is a simple 
geometric factor, depending on the choice of a. We were 
able to find, that essentially more exact expression (we 
call it as the "FSBE improved", or FSBEi approach) can 
be obtained with a — 2, i.e. when 



1 



-I 



df 



(30) 



Figs. [3] and [4] demonstrate that FSBEi turns out to be 
even more accurate and stable than FSBE in our sim- 
plified model example calculations. Unfortunately we 
were not able to obtain analytical derivatives dRi/drj 
for the radii from Eq. ( 30 ) for the specific surface imple- 
mentation we use. Nevertheless, the FSBE in the form 
presented here gives accurate enough for practical appli- 
cations values of solvation energies. Moreover in typical 
situations such as in proteins ions normally sit next to 
the water interfaces, and therefore, the resulting error 
for solvation energy is small. 

The idea to use integrals of the form of Eq. (29) in 



1 



[R(r) 



Jr w |r' - r 



Q! + 2 



df, 



(29) 



either volume or surface integral formulation to improve 
the accuracy of GB is not new [13j[27]. It was suggested 
that a linear combination of properly chosen integrals of 
the form of Eq. (29) with adjustable coefficients leads 
to a transferable (from small to big molecules) method. 
Nevertheless, such an approach does not let one to select 
a specific model (most of the models studied by the au- 
thors have similar errors when compared with the exact 
solution). We argue that FSBE method presented here 
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gives a unique approximation as a unique solution of the 
variational problem. 

FSBE has even more of subtle advantages over cur- 
rent GB approximations. We do not have exponential 
extrapolation factors in the denominator of Eq.(10) and 



thus are able to compute FSBE solvation energies consid- 
erably faster. FSBE lets us compute polarization surface 
charge density from Eq.(24) and hence obtain the sol- 



vation energy in essentially O(N) time and memory, as 
described in our subsequent work [10J. The deficiencies 
of the method, such as its (relative) failure to get large 
values of Born radii right, as well as its possible improve- 
ments, such as FSBEi, are left for future work. 

With all the apparent success of the method in 
solving the electrostatics problem, its applications to 
biomolecules modeling is limited by the fact, that wa- 
ter is not a simple dielectric with local and large value 
of the dielectric constant. The Poisson equation can de- 
scribe neither volume or surface phase transitions and 



hydrogen bonds networks rearrangements [8j |2T] nor wa- 
ter molecule orientational interactions in a polar liquid 
[9]. Nevertheless, the idea to prescribe Born approxima- 
tion a variational interpretation may serve as a universal 
framework to generate approximate solutions of arbitrary 
partial differential equations, including those of more so- 
phisticated water models, such as [7]. 
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